function M = extraRTS2array(filenin,delimiter,str_cell,N)
%extraRTS2array - extraction of topographical data provided by a total station  
%
%           M = extraRTS2array(FILENIN,DELIMITER,STR_CELL,N)
%
% This function extracts the topographical data from the ASCII file 
% FILENIN, assuming that they are in the form
%   [Date Reflector Angle1 Angle2 Distance x y z],
% where Date is the acquisition date in the form dd/mm/yy hh:mm, 
% Reflector is a string consistent with the input variable STR_CELL,
% Angle1, Angle2 and Distance are the polar coordinates and x, y, z
% are the Cartesian coordinates. If FILENIN is undefined or empty, the 
% filename can be interactively managed.
% DELIMITER is the delimiter of the ASCII file. If it is undefined or 
% empty, DELIMITER = '\t' (tabulation) is used. 
% STR_CELL is a cell variable of strings whith the target names.
% If STR_CELL is undefined or empty, this cell is used:
%   STR_CELL = {'P1','P2','P3','P4','P5','P6','P7','P8','P9','P10','P11',...
%     'P12','P13','P14','P15','P16','P17','P18','P19','P20','P21','P22',...
%     'P23','P24','P25','P26','P27','P28','P29','P30','P31',...
%     'N1','N2','N3','N4','N5','N6','N7','N8','N9','N10',...
%     'GPS_ROV 1','GPS_ROV 2',...
%     'R1','R1_BIS','R2','R2BIS','R3','R4','R5'}.
% The number N is the number of last considered rows. Is undefined or
% empty, it is N = 36000.
% The data related to the k-th target are placed in the matrix M(:,:,k) of
% the output 3D array M.
% The sampling rate is assumed to be 24 cpd. For the targets characterized 
% by sampling rate lower than 24 cpd, the positions for which no data are 
% provided are NaN.
%
% See also velArray, selectMeanVelocity.
%
%           M = extraRTS2array(FILENIN,DELIMITER,STR_CELL,N)

% G. Teza, 2020.

if nargin < 4
    N = 35000; 
end

if nargin < 3 || isempty(str_cell)
    str_cell = {'P1','P2','P3','P4','P5','P6','P7','P8','P9','P10','P11',...
        'P12','P13','P14','P15','P16','P17','P18','P19','P20','P21','P22',...
        'P23','P24','P25','P26','P27','P28','P29','P30','P31',...
        'N1','N2','N3','N4','N5','N6','N7','N8','N9','N10',...
        'GPS_ROV 1','GPS_ROV 2',...
        'R1','R1_BIS','R2','R2BIS','R3','R4','R5'};
end

if nargin < 2 || isempty(delimiter)
    delimiter = '\t';
end

if nargin < 1
    filenin = []; 
end
if isempty(filenin)
    [fil1,path1] = uigetfile(...
        {'*.txt','ASCII file (*.txt)'},...
        'INPUT FILE');
    if isnumeric(fil1) || isnumeric(path1)
        error('Invalid choice');
    else  
        filenin = fullfile(path1,fil1);
    end
end

FID = fopen(filenin);
C = textscan(FID,'%s%s%f%f%f%f%f%f','delimiter',delimiter);
fclose(FID);

c1 = C{1};
c2 = C{2};
c3 = C{3};
c4 = C{4};
c5 = C{5};
c6 = C{6};
c7 = C{7};
c8 = C{8};

% numerical data:
a = [c3 c4 c5 c6 c7 c8];

% to remove the useless rows
nr = size(c1,1);
Ie = zeros(nr,1); 
for k = 1:nr
    Ie(k) = isempty(c1{k});
end
Ie = logical(Ie);

c1(Ie) = [];
c2(Ie) = [];
a(Ie,:) = [];

c1 = c1(end-N+1:end,:);
c2 = c2(end-N+1:end,:);
a = a(end-N+1:end,:);

nrow = size(a,1);

ns = numel(str_cell);   % number of stations

c11 = c1(1);
d1 = obsdate(c11);              % inital date
c1end = c1(end);
dend = obsdate(c1end);        % final date

dv = (floor(d1):(1/24):ceil(dend))';    % date vector
nm = length(dv);
M  = nan(nm,7,ns);      % output matrix inizialization
for m = 1:ns
    M(:,1,m) = dv;
end

dateb = zeros(nrow,1);
for h = 1:nrow
    dateb(h) = obsdate(c1(h));
    datestr(dateb(h));
end

for k = 1:ns
    strk = str_cell{k};
    fprintf('Processing of station %s in progress\n',strk);
    Ik = strcmpi(strtrim(c2),strk); % search of rows related to k-th station
    amk = a(Ik,:);                  % related data matrix
    bvk = dateb(Ik,1);              % related time vector
    nk = length(bvk);
    for h = 1:nk
        th = bvk(h);
        [minh,Ih] = min(abs(dv-th));   % search of the observation time in dh
        if minh < 1/24              % to avoid data of other acquisition times 
            M(Ih,2:7,k) = amk(h,:);
        end
    end
end

function dk = obsdate(bk)
bk = char(bk);
dd = str2double(bk(1:2));
mm = str2double(bk(4:5));
yy = str2double(bk(7:10));
if length(bk) > 10
    h = str2double(bk(12:13));
    m = str2double(bk(15:16));
else
    h = 0;
    m = 0;
end
dk = datenum(yy,mm,dd,h,m,0);